El objetivo de este informe es presentar un protocolo en R, utilizando el paquete clusterProfiler de Bioconductor, para realizar el análisis de significación biológica de una lista de genes. Este protocolo se utilizará como referencia para el desarrollo de la aplicación Shiny objeto del TFM.
El análisis de significación biologica permite identificar aquellos procesos biológicos, como por ejemplo características funcionales, vías metabólicas o cascadas de transducción de señales, que están relacionados con los datos contenidos de la lista de genes de partida. El análisis se puede dividir en tres etapas principales:
Las anotaciones contenidas en bases de datos como GO, KEGG o Reactome permiten relacionar conjuntos de genes que trabajan coordinadamente en determinados procesos biológicos, asignándolos a ciertas categorías o términos. A pesar de que existen numerosas bases de datos para realizar este proceso, las que se acaban de citar suelen ser las primeras opciones a la hora de realizar estos análisis debido a su larga tradición, actualización y diversidad de especies disponibles.
GO: El proyecto Gene Ontology es una iniciativa bioinformática que tiene como objetivo estandarizar la representación de los genes y los atributos de sus productos génicos de todas las especies. Los términos indicados están relacionados unos con otros de manera jerárquica, desde términos más generales hacia términos más específicos. Las anotaciones incluyen información sobre Proceso biológico (BP): eventos celulares a los que contribuye el producto génico, Función molecular (MF): descripción bioquímica del producto génico y Componente celular (CC): localización o complejos de los que forma parte el gen o su producto génico.
KEGG: La Kyoto Encyclopedia of Genes and Genomes es un conjunto de base de datos diseñadas para facilitar la comprensión de sistemas biológicos (células, organismos o ecosistemas) a nivel molecular. Su base de datos más conocida es KEGG PATHWAY, una colección de mapas de procesos bioquímicos que representan interacciones moleculares, reacciones metabólicas, procesos celulares y relacionados con enfermedades, entre otros.
Reactome: Reactome es una base de datos de reacciones, vías y procesos biológicos, disponible en línea y de código abierto. El modelo de datos de Reactome considera las reacciones entre ácidos nucleicos, proteínas, complejos y moléculas pequeñas, formando una red de interacciones biológicas que se agrupan en categorías. Los ejemplos de rutas biológicas en Reactome incluyen vías de señalización, regulación transcripcional, traducción, apoptosis y metabolismo, entre otras.
Una vez hemos asignado los genes de nuestra lista a ciertas
categorías, debemos analizar cuáles de estas anotaciones están
relacionadas con el problema que se está estudiando o bien aparecen por
azar entre las muchas anotaciones de los genes de la lista. Actualmente
existen numerosos métodos destinados a realizar este análisis. En
general, los métodos más utilizados se basan en el análisis de
enriquecimiento, en el que se establece si una determinada categoría
aparece más o menos a menudo en la lista de genes seleccionados respecto
a una población (universo génico) de referencia. De manera simple, se
considera que si una categoría aparece más a menudo en nuestra lista que
en la de referencia es probable que sea relevante para explicar las
diferencias observadas.
La nomenclatura y clasificación de los métodos de análisis propuesto
por los expertos en la materia es diversa. De acuerdo a la utilizada por
Khatri et al., se agrupan los métodos de análisis en diferentes
generaciones, distinguiendo los de primera generación (ORA: Over
Representation Analysis), segunda generación (FCS: Functional
Class Analysis) y tercera generación (PT: Pathway
Topology). Los métodos ORA y FCS consideran únicamente el número de
genes correspondientes a un proceso biológico determinado, mientras que
los métodos PT tienen en cuenta la información adicional relacionada con
cómo y dónde los genes pueden interaccionar entre ellos. Relacionado con
esta distinción, existen otras clasificaciones propuestas, como la de
Geistlinger en la que se distinguen dos grandes grupos de
métodos basados en el análisis de enriquecimiento: Set-based
enrichment analysis, en el que no se considera la posible
interacción entre genes, y el Network-based enrichment
analysis, en el que se consideran las anotaciones relacionadas con
estas interacciones.
A pesar de que los métodos de última generación (Pathway Topology o Network-based enrichment analysis) parecen aproximar mejor la realidad de los procesos biológicos, existen numerosas limitaciones a nivel de anotación y metodología. Por ello, el presente TFM se centrará en los dos métodos de enriquecimiento más utilizados:
ORA (Over Representation Analysis) o Análisis de sobre-representación: Comprueba la sobre-representación estadística de una lista de genes de interés en una lista de referencia.
GSEA (Gene Set Enrichment Analysis) o Análisis de enriquecimiento de conjuntos de genes: Incorpora los valores de expresión y estadísticos diferenciales de todos los genes de para hacer un test de análisis de significación estadística, con el objetivo de comprobar si los genes correspondientes a ciertas categorías se acumulan en la parte superior o inferior de la lista ordenada por dirección y magnitud del cambio de expresión.
Durante el análisis del presente trabajo, se utilizan los métodos anteriores anotando los genes con Go, KEGG y Reactome. Por tanto, en este informe se presentan seis análisis diferenciados: ORA-GO, ORA-KEGG, ORA-Reactome, GSEA-GO, GSEA-KEGG y GSEA-Reactome.
Por último, se evaluará el uso de un método de tercera generación:
La visualización de los resultados del análisis permite identificar e interpretar los procesos biológicos enriquecidos en nuestro análisis. Además, muchas de las técnicas de visualización utilizadas permiten minimizar redundancias que frecuentemente aparecen en los resultados de enriquecimiento. La visualización permitirá agrupar procesos y vías similares en categorías funcionales comunes. Los paquetes clusterProfiler, enrichplot, Pathview, ReactomePA y SPIA contienen algunos de los métodos de visualización más utilizados, entre los que destacan los siguientes:
Barplot: La función barplot() muestra un gráfico de barras para visualizar las categorías enriquecidas. El color representa estadísticos de enriquecimiento (por ejemplo, valores p) y la longitud de la barra respecto el eje de abcisas el número de genes, o proporción respecto al total, anotados en la categoría.
Dotplot: La función dotplot() es similar a la barplot(), pero puede representar una característica más por el tamaño del punto. El color se relaciona con los valores p ajustados del análisis y el tamaño del punto con el número de genes anotados a la categoría. En el eje de abcisas número de genes, o proporción respecto al total, anotados en la categoría.
Enrichment Map: El mapa de enriquecimiento organiza las categorías enriquecidas en una red donde se conectan conjuntos de genes solapados. Los conjuntos de genes que se solapan tienden a agruparse, lo que facilita la identificación de módulos funcionales. Para visualizar mapas de enriquecimiento tenemos diferentes opciones: la función emmaplot(), aplicable en todas las anotaciones utilizadas y la función goplot(), específica para anotaciones en GO.
Category Netplot: La función cnetplot() representa los vínculos entre genes y categorías biológicas (por ejemplo, términos GO o vías KEGG) en forma de red. El resultado de GSEA también es compatible pero solo se muestran los genes enriquecidos en el núcleo.
Upsetplot: La función upsetplot() es una alternativa al Category Netplot para visualizar la asociación compleja entre genes y conjuntos de genes. Esta función destaca el solapamiento de genes entre diferentes categorías.
Ridgeline plot: La función ridgeplot() visualiza las distribuciones de expresión para las categorías obtenidas en el análisis GSEA. Ayuda a los usuarios a interpretar las vías biológicas enriquecidas que derivan de sobreexpresión o infraexpresión de los conjuntos de genes identificados.
GSEA plot: La función gseaplot() (también gseaplot2()) del paquete enrichplot permite representar para el conjunto de genes seleccionado el Running Enrichment score en la lista ordenada de acuerdo al análisis GSEA.
Pathview: La función pathview() del paquete Pathview es una herramienta que permite integrar y visualizar datos de vías biológicas. Los usuarios solo necesitan especificar la vía deseada, en base al análisis de enriquecimiento, y se representa un gráfico con las vistas nativas de KEGG.
viewPath: La función viewPath() del paquete permite representar la red que conecta los genes relacionados con una categoría especificada de Reactome.
plotP: La función plotP genera un gráfico de significación estadística de las vías de señalización biológica obtenidas en el análisis realizado con SPIA. Las líneas oblicuas del gráfico muestran los umbrales correspondientes a regiones de significancia basada en diferentes criterios. Los puntos representados a la derecha de estos umbrales corresponden a catagerías KEGG significativamente enriquecidas.
Se ha creado un proyecto de RStudio en la carpeta principal de
trabajo. Se ha creado la carpeta Data donde se ubicará la lista
de genes utilizada como punto de partida del análisis y la carpeta
Results para guardar los archivos generados.
Para ejecutar el código se necesita tener instalado Bioconductor.
> if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager")
Además, se necesitarán instalar los siguientes paquetes:
> if (!require(clusterProfiler)) BiocManager::install("clusterProfiler")
> if (!require(pathview)) BiocManager::install("pathview")
> if (!require(enrichplot)) BiocManager::install("enrichplot")
> if (!require(ReactomePA)) BiocManager::install("ReactomePA")
> if (!require(DOSE)) BiocManager::install("DOSE")
> if (!require(SPIA)) BiocManager::install("SPIA")
> if (!require(ggplot2)) install.packages("ggplot2", dep=TRUE)
> if (!require(ggridges)) install.packages("ggridges", dep=TRUE)
> if (!require(ggupset)) install.packages("ggupset", dep=TRUE)
En base al organismo origen de los datos de los genes en evaluación,
se deberá descargar el paquete de anotaciones correspondiente. El
usuario deberá seleccionar el organismo, teniendo en cuenta el código
requerido según la ontología o base de datos utilizada.
GO: Considerar los paquestes actualizados a utilizar en el siguiente
enlace: OrgDb.
KEGG: Considerar la codificación a utilizar en el siguiente enlace: KEGG.
Reactome: Considerar la codificación a utilizar en el siguiente enlace:
Reactome
A continuación definimos unas variables relacionadas con el organismo en estudio según la base de datos de anotaciones.
> #Seleccionamos el organismo relacionado con los datos de partida. En este caso Ratón.
> #Anotación en GO
> organismo_go = "org.Mm.eg.db" #En la aplicación se permitirá seleccionar Humano, ratón y Rata
> BiocManager::install(organismo_go, character.only = TRUE)
> library(organismo_go, character.only = TRUE)
> #Anotación en KEGG
> organismo_kegg = "mmu"
> #Anotación en Reactome
> organismo_reactome="mouse"
Por último, cargamos los paquetes requeridos para el
análisis.
> library(clusterProfiler)
> library(enrichplot)
> library(ReactomePA)
> library(ggplot2)
> library(DOSE)
> library(pathview)
> library(SPIA)
Los datos de partida de este trabajo son una lista de genes diferencialmente expresados (topTable), obtenidos a partir de un estudio de comparación de la expresión génica en células luminales y basales extraídas de glándulas mamarias de ratones vírgenes, gestantes de 18,5 días y lactantes de 2 días. Los datos del estudio están disponibles en Gene Expression Omnibus (GEO) con el número de acceso GSE60450.
La lista de genes sobre la que se realiza el análisis de significación biológica es la selección de genes obtenida tras realizar el análisis de expresión diferencial con limma-voom entre los grupos gestantes y lactantes: toptab_B.PregVsLac.csv. Obtenida del análisis de expresión diferencial relizado por Mireia Ferrer et al. Esta lista contiene los genes identificados con el identificador ENTREZ en la primera columna y varios parámetros estadísticos obtenidos en el análisis de expresión diferencia para cada gen, relacionados con la diferencia entre los estados “pregnant” y “lactate” en el grupo “Basal”:
Antes de realizar el análisis, el protocolo requiere la estandarización de la siguiente información de la lista de genes:
Opcionalmente, la aplicación Shiny en desarrollo permitirá transformar los datos de partida de manera automática para cumplir los requerimientos anteriores. También se permitirán cargar los datos de archivos .csv, .txt o .xlsx.
Cargamos la lista de genes y preparamos los datos.
> df = read.csv("Data/toptab_B.PregVsLac.csv", header=TRUE)
> str(df) #Comprobamos que la tabla es un data frame
'data.frame': 15804 obs. of 7 variables:
$ X : int 24117 381290 78896 226101 16012 231830 16669 55987 231991 14620 ...
$ logFC : num 1.82 -2.14 2.81 -2.33 -2.9 ...
$ AveExpr : num 2.98 3.94 3.04 6.22 1.98 ...
$ t : num 20.1 -19.1 18.5 -18.3 -18.2 ...
$ P.Value : num 1.06e-10 1.98e-10 2.76e-10 3.30e-10 3.41e-10 ...
$ adj.P.Val: num 1.02e-06 1.02e-06 1.02e-06 1.02e-06 1.02e-06 ...
$ B : num 15 14.4 14.1 13.9 13.5 ...
> head(df) #Comprobamos los nombres de las diferentes variables
X logFC AveExpr t P.Value adj.P.Val B
1 24117 1.819943 2.975545 20.10780 1.063770e-10 1.01624e-06 14.96977
2 381290 -2.143885 3.944066 -19.07495 1.982934e-10 1.01624e-06 14.39556
3 78896 2.807548 3.036519 18.54773 2.758828e-10 1.01624e-06 14.07416
4 226101 -2.329744 6.223525 -18.26861 3.297667e-10 1.01624e-06 13.85802
5 16012 -2.896115 1.978449 -18.21525 3.413066e-10 1.01624e-06 13.46984
6 231830 2.253400 4.760597 18.02627 3.858161e-10 1.01624e-06 13.67600
> #En la aplicación se mostrarán los requerimientos de los datos de partida a nivel de identificadores de genes
> #y nombres de los parámetros estadísticos a evaluar.
> #En caso necesario, se realizará la transformación de los identificadores con la función bitr()
> #bitr(geneID, fromType, toType, OrgDb, drop = TRUE)
> #La primera columna debe corresponder a los identificadores de los genes.
> colnames(df)[1]<-"EntrezID" #Llamamos EntrezID a la variable que contiene los EntrezID
El análisis de sobrerrepresentación (ORA) es un método estadístico que determina si los genes correspondientes a una categoría determinada están sobrerrepresentados en un subconjunto de nuestros datos. Para determinar si alguna categoría está sobrerrepresentada, se determina la probabilidad de tener la proporción observada de genes asociados a esta categoría en nuestro subconjunto respecto a la proporción de genes asociados con la misma categoría en la lista de referencia (“universo”).
Para evaluar si alguna categoría está sobrerrepresentada en el subconjunto, se puede utilizar la distribución hipergeométrica. Esto corresponde a realizar el test exacto de Fischer de una cola. A partir de este análisis se obtienen unos valores p estadísticos para cada categoría que deberán ajustarse por un método de comparaciones múltiples (p.adjust) y para el control de la tasa de falsos positivos FDR (valor q).
En este caso, el subconjunto en el que evaluaremos la sobrerrepresentación es una lista de genes obtenida a partir de de la original en la que hemos aplicado unos filtros relacionados con los estadísticos y medidas de cambio de expresión génica obtenidos en el estudio de expresión diferencial. En el ejemplo, el subconjunto está formado por aquellos genes con un adj.P.val <0.05 y logFC>2 (up-regulated). En la aplicación en desarrollo, el usuario podrá seleccionar el umbral deseado para el análisis, permitiendo realizar el análisis de los genes sobre o subexpresados por separado, o ambos a la vez (abs(logFC)>2)).
Opcionalmente, la aplicación permitirá seleccionar como la lista de genes a evaluar la lista inicial sin filtrar. Dependiendo de la necesidad del análisis, el universo de genes utilizado en el análisis podrá ser la misma topTable inicial o bien todo el genoma del organismo.
Preparamos los datos sobre los que realizaremos el análisis.
> # Generamos la lista de genes de referencia (universo) con la que se comparará nuestra lista
> # En este caso, el universo será la lista de genes original (topTable inicial)
> universe_list <- as.character(df$EntrezID)
> # Eliminamos posibles valores NA
> universe_list<-na.omit(universe_list)
> # Eliminamos posibles valores duplicados
> universe_list <- universe_list[which(duplicated(universe_list) == F)]
>
> # Otra opción sería utilizar como universo todos los genes del genoma del organismo
> # En ese caso, se requería la siguiente instrución:
> # universe_all_list <- select(org.Mm.eg.db, keys=keys(org.Mm.eg.db), columns="ENTREZID")
> # universe_list <- as.character(universe_all_list$ENTREZID)
>
> # De la lista de genes original, seleccionamos los de interés en base a valor adj.P.val (padj < 0.05)
> # En la aplicación se añadirán selectores para filtrar los datos por adj.P.Val y logFC
> sig_genes_df = subset(df, adj.P.Val < 0.05)
> # De la lista anterior, extraemos los valores de logFC para filtrar posteriormente
> ora_list <- sig_genes_df$logFC
> # Nombramos el vector anterior con los identificadores
> names(ora_list) <- sig_genes_df$EntrezID
> # omitimos NA values
> ora_list <- na.omit(ora_list)
> # Eliminamos posibles valores duplicados
> ora_list <- ora_list[which(duplicated(names(ora_list)) == F)]
> # Finalmente filtramos por el valor de logFC deseado. En este caso nos centraremos en up-regulated (logFC > 2)
> ora_list <- ora_list[ora_list > 2]
En el análisis ORA se pueden personalizar algunos parámetros independientemente de la base de datos de anotaciones utilizada. Los más destacados son:
minGSSize: Número mínimo de genes para cada conjunto de genes (categoría). Si es inferior al establecido, esa categoría no se reporta en los resultados.
maxGSize: Nmero máximo de genes para cada conjunto de genes (categoría). Si es superior al establecido, esa categoría no se reporta en los resultados.
ont: Ontología GO. Opciones: “BP”, “MF”, “CC” o “ALL”.
pAdjustMethod: Método de ajuste para comparaciones múltiples: “holm”, “hochberg”, “hommel”, “bonferroni”, “BH”, “BY”, “fdr”, “none”.
pvalueCutoff: p-valor de corte. Sólo se reportan las categorías con un p-valor inferior al establecido.
qvalueCutoff: q-valor de corte. Sólo se reportan las categorías con un q-valor inferior al establecido.
Realizamos el análisis con la función enrichGO(). Se debe indicar en el argumento gene el subconjunto de genes en evaluación y en argumento universe la lista de genes de referencia. Esta función admite numerosos gene ID. Se pueden consultar con la siguiente función (keytypes()). Es importante confirmar que el identificador utilizado en nuestras listas de genes esté soportado por la función y la anotación seleccionada. En nuestro caso seguimos trabajando con ENTREZID.
> keytypes(org.Mm.eg.db)
[1] "ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT" "ENSEMBLTRANS"
[6] "ENTREZID" "ENZYME" "EVIDENCE" "EVIDENCEALL" "GENENAME"
[11] "GENETYPE" "GO" "GOALL" "IPI" "MGI"
[16] "ONTOLOGY" "ONTOLOGYALL" "PATH" "PFAM" "PMID"
[21] "PROSITE" "REFSEQ" "SYMBOL" "UNIPROT"
El objeto resultante almacena las categorías GO enriquecidas en el subconjunto de genes, la proporción observada en el subconjunto y en la lista de referencia, y los valores estadísticos que justifican que dichas categorías se encuentran sobrerrepresentadas.
> #Seleccionamos la ontología BP de GO.
> ora_go <- enrichGO(gene = names(ora_list),
+ universe = universe_list,
+ OrgDb = organismo_go,
+ keyType = "ENTREZID", #Se podría trabajar también con ID Ensemble
+ readable = T,
+ ont = "BP",
+ pAdjustMethod = "BH",
+ qvalueCutoff = 0.05)
Como se ha comentado anteriormente, las categorías anotadas en GO están relacionada de manera jerárquica, por lo que un los procesos relacionados con un término padre puede superponerse en gran medida con los de sus hijos. Esto puede producir resultados redundantes, en los que tengamos muchas categorías relacionadas con el mismo proceso biológico. Para solucionar este problema, clusterProfiler implementa un método de simplificación con la función simplify() para reducir las categorías redundantes de GO obtenidas tanto en el análisis ORA como GSEA. Existen otros paquetes que ofrecen funcionalidades similares, e incluso con opciones de visualización avanzadas de las categorías solapadas, como simplifyEnrichment o rrvgo. No obstante, con el objetivo de optimizar recursos en la aplicación y simplificar el análisis, utilizaremos simplify() del paquete clusterProfiler.
En los argumentos del método simplify() se aplica select_fun (que puede ser una función definida por el usuario) a la característica by para seleccionar un término representativo de entre los términos redundantes (que tienen una similitud mayor que el umbral de corte, cutoff).
> ora_go_simplify <- simplify(ora_go, cutoff=0.7, by="p.adjust", select_fun=min)
En las siguientes instrucciones utilizaremos el objeto original del análisis (ora_go) a modo de ejemplo, pero se podría trabajar indistintamente con el original o el simplificado en base al criterio del usuario.
Podemos guardar los resultados obtenidos del análisis ORA-GO en un archivo csv.
> ora_go_results <- data.frame(ora_go)
> head(ora_go_results)[,2:7]
Description GeneRatio BgRatio
GO:0007059 chromosome segregation 28/282 326/14567
GO:0000070 mitotic sister chromatid segregation 21/282 180/14567
GO:0000819 sister chromatid segregation 21/282 204/14567
GO:0140014 mitotic nuclear division 24/282 288/14567
GO:0098813 nuclear chromosome segregation 23/282 266/14567
GO:0051276 chromosome organization 31/282 475/14567
pvalue p.adjust qvalue
GO:0007059 3.921301e-11 7.418301e-08 6.171278e-08
GO:0000070 4.256054e-11 7.418301e-08 6.171278e-08
GO:0000819 4.507580e-10 5.237808e-07 4.357328e-07
GO:0140014 1.914910e-09 1.444053e-06 1.201306e-06
GO:0098813 2.071217e-09 1.444053e-06 1.201306e-06
GO:0051276 3.194124e-09 1.855786e-06 1.543827e-06
> write.csv(ora_go_results, "Results/clusterProfiler_ORA_GO.csv")
Visualización del Upset plot:
> upsetplot(ora_go)
Visualización del Barplot:
> barplot(ora_go,
+ drop = TRUE,
+ showCategory = 10,
+ title = "GO Biological Pathways",
+ font.size = 9)
Visualización del Dotplot:
> dotplot(ora_go, showCategory = 10, font.size=9)
Visualización del Enrichment map:
> ora_go_sim <- pairwise_termsim(ora_go)
> emapplot(ora_go_sim, showCategory=10, color="pvalue", cluster.params=list(cluster=TRUE, legend=TRUE))
Visualización del GO plot:
> goplot(ora_go, showCategory = 5, cex=0.5)
Visualización del Category Netplot
> cnetplot(ora_go, showCategory=15, color.params =list(foldChange =ora_list))
Para realizar el análisis ORA con anotación de las vías KEGG utilizaremos la función enrichKEGG(). En este caso es importante confirmar que nuestra lista de genes utiliza el identificador ENTREZID, ya que es uno de los identificadores compatibles con la función enrichKEGG(). En caso de que nuestra lista de genes no esté identificada con ENTREZID, se tendrá que realizar la conversión utilizando la función bitr() de clusterProfiler.
Como se ha comentado anteriormente, en este trabajo la lista de genes de partida está codificada en ENTREZ ID, por lo que no se requiere de ninguna acción relacionada.
El objeto resultante almacena las vías KEGG enriquecidas en el subconjunto de genes, la proporción observada en el subconjunto y en la lista de referencia, y los valores estadísticos que justifican que dichas categorías se encuentran sobrerrepresentadas.
> ora_kegg <- enrichKEGG(gene=names(ora_list),
+ universe=universe_list,
+ organism=organismo_kegg,
+ pAdjustMethod = "BH",
+ qvalueCutoff = 0.05,
+ keyType = "ncbi-geneid")
Podemos guardar los resultados obtenidos del análisis ORA-KEGG en un archivo csv.
> ora_kegg_results <- data.frame(ora_kegg)
> head(ora_kegg_results)[2:7]
Description
mmu05150 Staphylococcus aureus infection - Mus musculus (house mouse)
mmu04110 Cell cycle - Mus musculus (house mouse)
mmu04080 Neuroactive ligand-receptor interaction - Mus musculus (house mouse)
mmu05322 Systemic lupus erythematosus - Mus musculus (house mouse)
mmu05323 Rheumatoid arthritis - Mus musculus (house mouse)
mmu04670 Leukocyte transendothelial migration - Mus musculus (house mouse)
GeneRatio BgRatio pvalue p.adjust qvalue
mmu05150 10/138 47/6148 5.975321e-08 1.434077e-05 1.264252e-05
mmu04110 14/138 151/6148 6.080825e-06 7.296990e-04 6.432873e-04
mmu04080 12/138 145/6148 8.928915e-05 4.838974e-03 4.265937e-03
mmu05322 8/138 65/6148 9.027976e-05 4.838974e-03 4.265937e-03
mmu05323 8/138 66/6148 1.008120e-04 4.838974e-03 4.265937e-03
mmu04670 9/138 91/6148 1.833358e-04 7.333432e-03 6.464999e-03
> write.csv(ora_kegg_results, "Results/clusterProfiler_ORA_KEGG.csv")
Visualización del Upset plot:
> upsetplot(ora_kegg)
Visualización del Barplot:
> barplot(ora_kegg,
+ showCategory = 10,
+ title = "Enriched Pathways",
+ font.size = 8)
Visualización del Dotplot:
> dotplot(ora_kegg,
+ showCategory = 10,
+ title = "Enriched Pathways",
+ font.size = 8)
Visualización del Enrichment map:
> ora_kegg_sim <- pairwise_termsim(ora_kegg)
> emapplot(ora_kegg_sim, showCategory=10, color="pvalue", cluster.params=list(cluster=TRUE, legend=TRUE))
Visualización del Category Netplot:
> cnetplot(ora_kegg, showCategory=5, color.params =list(foldChange =ora_list))
Visualización del Pathview: Estas instrucciones crearán diferentes archivos (PNG y PDF) con información de las vías KEGG enriquecidas. El usuario deberá seleccionar qué vía de las obtenidas en el análisis quiere visualizar.
> head(ora_kegg)[,1:2]
ID
mmu05150 mmu05150
mmu04110 mmu04110
mmu04080 mmu04080
mmu05322 mmu05322
mmu05323 mmu05323
mmu04670 mmu04670
Description
mmu05150 Staphylococcus aureus infection - Mus musculus (house mouse)
mmu04110 Cell cycle - Mus musculus (house mouse)
mmu04080 Neuroactive ligand-receptor interaction - Mus musculus (house mouse)
mmu05322 Systemic lupus erythematosus - Mus musculus (house mouse)
mmu05323 Rheumatoid arthritis - Mus musculus (house mouse)
mmu04670 Leukocyte transendothelial migration - Mus musculus (house mouse)
> dme <- pathview(gene.data=universe_list, pathway.id=ora_kegg@result$ID[1], species = organismo_kegg)
También se puede navegar directamente a la página web de KEGG para visualizar la vía del término seleccionado.
> browseKEGG(ora_kegg, ora_kegg@result$ID[1])
ORA-KEGG Pathway para mmu05150
Para realizar el análisis ORA con anotación con los términos Reactome utilizaremos la función enrichPathway() del paquete ReactomePA. De nuevo, la función requiere que nuestra lista de genes sea un vector de identificadores ENTREZ. En caso de que nuestra lista de genes no esté identificada con ENTREZ ID, se tendrá que realizar la conversión utilizando la función bitr() de clusterProfiler.
Como se ha comentado anteriormente, en este trabajo la lista de genes de partida está codificada en ENTREZ ID, por lo que no se requiere de ninguna acción relacionada.
El objeto resultante almacena las vías Reactome enriquecidas en el subconjunto de genes, la proporción observada en el subconjunto y en la lista de referencia, y los valores estadísticos que justifican que dichas categorías se encuentran sobrerrepresentadas.
> ora_reactome <- enrichPathway(gene=names(ora_list),
+ universe=universe_list,
+ organism=organismo_reactome,
+ pAdjustMethod = "BH",
+ pvalueCutoff = 0.05,
+ readable = TRUE)
Podemos guardar los resultados obtenidos del análisis ORA-REACTOME en un archivo csv.
> ora_reactome_results <- data.frame(ora_reactome)
> head(ora_reactome_results)[2:7]
Description
R-MMU-2500257 Resolution of Sister Chromatid Cohesion
R-MMU-141424 Amplification of signal from the kinetochores
R-MMU-141444 Amplification of signal from unattached kinetochores via a MAD2 inhibitory signal
R-MMU-69618 Mitotic Spindle Checkpoint
R-MMU-69620 Cell Cycle Checkpoints
R-MMU-9648025 EML4 and NUDC in mitotic spindle formation
GeneRatio BgRatio pvalue p.adjust qvalue
R-MMU-2500257 16/153 111/6954 1.858618e-09 7.341542e-07 6.417124e-07
R-MMU-141424 13/153 88/6954 4.920479e-08 5.542674e-06 4.844762e-06
R-MMU-141444 13/153 88/6954 4.920479e-08 5.542674e-06 4.844762e-06
R-MMU-69618 14/153 105/6954 5.612834e-08 5.542674e-06 4.844762e-06
R-MMU-69620 21/153 257/6954 1.639365e-07 1.295098e-05 1.132025e-05
R-MMU-9648025 13/153 100/6954 2.321470e-07 1.528301e-05 1.335863e-05
> write.csv(ora_reactome_results, "Results/clusterProfiler_ORA_reactome.csv")
Visualización del Upset plot:
> upsetplot(ora_reactome)
Visualización del Barplot:
> barplot(ora_reactome,
+ showCategory = 10,
+ title = "Enriched Pathways",
+ font.size = 8)
Visualización del Dotplot:
> dotplot(ora_reactome,
+ showCategory = 10,
+ title = "Enriched Pathways",
+ font.size = 8)
Visualización del Enrichment map:
> ora_reactome_sim <- pairwise_termsim(ora_reactome)
> emapplot(ora_reactome_sim, showCategory=10, color="pvalue", cluster.params=list(cluster=TRUE, legend=TRUE))
Visualización del Category Netplot:
> cnetplot(ora_reactome, showCategory=5, color.params =list(foldChange =ora_list))
Visualización del viewPathway: Esta función permite visualizar la via Reactome indicada por el usuario.
> head(ora_reactome)[,1:2]
ID
R-MMU-2500257 R-MMU-2500257
R-MMU-141424 R-MMU-141424
R-MMU-141444 R-MMU-141444
R-MMU-69618 R-MMU-69618
R-MMU-69620 R-MMU-69620
R-MMU-9648025 R-MMU-9648025
Description
R-MMU-2500257 Resolution of Sister Chromatid Cohesion
R-MMU-141424 Amplification of signal from the kinetochores
R-MMU-141444 Amplification of signal from unattached kinetochores via a MAD2 inhibitory signal
R-MMU-69618 Mitotic Spindle Checkpoint
R-MMU-69620 Cell Cycle Checkpoints
R-MMU-9648025 EML4 and NUDC in mitotic spindle formation
> viewPathway(ora_reactome@result$Description[1], readable=TRUE)
A diferencia del análisis ORA, el análisis de enriquecimiento de conjunto de genes (GSEA) utiliza todos los genes de la lista de partida. Normalmente, los genes se ordenan en base a valores estadísticos o medidas de diferencia de expresión entre los grupos en evaluación, en nuestro caso logFC, y a continuación, se analiza si las categorías anotadas se distribuyen aleatoriamente por toda la lista de genes o se encuentran principalmente en la parte superior o inferior.
Hay tres elementos claves en el análisis GSEA:
Cálculo de la puntuación de enriquecimiento (ES: enrichmentScore): Esta puntuación representa la magnitud en que una categoría está sobrerrepresentada en la parte superior o inferior de la lista ordenada de genes.
Estimación del nivel de significación del ES: El p-valor del ES se calcula mediante una prueba de permutación.
Ajuste para pruebas de comparaciones múltiples: El p-valor se corrige de acuerdo al ajuste de pruebas de comparaciones múltiples (p.adjust) y para el control de la tasa de falsos positivos (q-valor).
Preparamos los datos sobre los que realizaremos el análisis.
> # Seleccionamos los valores de logFC de toda la lista de genes
> original_gene_list <- df$logFC
>
> # Nombramos el vector anterior con los gene ID
> names(original_gene_list) <- df$EntrezID
>
> # Eliminamos posibles valores NA
> gsea_list<-na.omit(original_gene_list)
>
> # Eliminamos posibles valores duplicados
> gsea_list <- gsea_list[which(duplicated(names(gsea_list)) == F)]
>
> # Ordenamos la lista en orden decreciente de logFC (requerido en ClusterProfiler)
> gsea_list = sort(gsea_list, decreasing = TRUE)
En el análisis GSEA se pueden configurar una serie de parámetros independientemente de la base de datos de anotaciones utilizada:
nPerm: Número de permutaciones utilizado en el análisis de significancia de la puntuación ES de cada categoría. A mayor número de permutaciones, el resultado obtenido será más exacto pero el análisis requerirá más tiempo de computación. No obstante, por recomendación de la documentación relacionada con esta función, no se suele indicar el número de permutaciones.
minGSSize: Número mínimo de genes para cada conjunto (categoría). Si es inferior al establecido, esa categoría no se reporta en los resultados.
maxGSSize: Número máximo de genes para cada conjunto (categoría). Si es superior al establecido, esa categoría no se reporta en los resultados.
pAdjustMethod: Ajuste del p-valor por el método de comparaciones múltiples seleccionado: “holm”, “hochberg”, “hommel”, “bonferroni”, “BH”, “BY”, “fdr”, “none”.
pvalueCutoff: p-valor de corte. Sólo se reportan las categorías con un p-valor inferior al establecido.
eps: Este parámetro establece el límite inferior para calcular el p-valor. Por defecto es 1e-10, pero se recomienda indicar 0 porque algunas vías presentan p-valores inferiores.
Realizamos el análisis con la función gseGO(). Seleccionamos toda nuestra lista de genes y configuramos diferentes parámetros del análisis. La aplicación Shiny permitirá configurar alguno de estos parámetros.
El objeto resultante almacena las categorías GO enriquecidas, la cuenta de genes anotados en ellas y los valores estadísticos que justifican que dichas categorías se encuentran significativamente enriquecidas.
> #Seleccionamos la ontología BP de GO. Opciones: “BP”, “MF”, “CC” o “ALL”
> gsea_go <- gseGO(geneList=gsea_list,
+ ont ="BP",
+ keyType = "ENTREZID",
+ minGSSize = 20,
+ pvalueCutoff = 0.05,
+ verbose =FALSE,
+ eps=0,
+ OrgDb = organismo_go,
+ pAdjustMethod = "BH")
De la misma manera que se ha realizado para el análisis ORA, la función simplify() nos permite reducir las categorías redundantes de GO en el análisis GSEA.
> gsea_go_simplify <- simplify(gsea_go, cutoff=0.7, by="p.adjust", select_fun=min)
Podemos guardar los resultados obtenidos del análisis GSEA-GO en un archivo csv.
> gsea_go_results <- data.frame(gsea_go)
> head(gsea_go_results)[2:8]
Description setSize enrichmentScore
GO:0000070 mitotic sister chromatid segregation 180 0.6404028
GO:0000819 sister chromatid segregation 204 0.6168462
GO:0007059 chromosome segregation 326 0.5535793
GO:0042254 ribosome biogenesis 295 0.5589020
GO:0022613 ribonucleoprotein complex biogenesis 408 0.5082035
GO:0098813 nuclear chromosome segregation 266 0.5581239
NES pvalue p.adjust qvalue
GO:0000070 2.430429 1.637096e-18 4.419766e-15 3.583515e-15
GO:0000819 2.380431 3.325633e-18 4.419766e-15 3.583515e-15
GO:0007059 2.232800 3.229983e-18 4.419766e-15 3.583515e-15
GO:0042254 2.232364 1.783819e-17 1.778022e-14 1.441607e-14
GO:0022613 2.080452 3.836421e-16 3.059162e-13 2.480347e-13
GO:0098813 2.212462 1.403784e-15 9.328142e-13 7.563192e-13
> write.csv(gsea_go_results, "Results/clusterProfiler_GSEA_GO.csv")
Visualización del Upset plot:
> upsetplot(gsea_go)
Visualización del Dotplot.
> dotplot(gsea_go, showCategory=10, font.size=8)
> # También podemos separar las categorías que provienen de sobre/subexpresión
> dotplot(gsea_go, showCategory=8, font.size=8, split=".sign") + facet_grid(.~.sign)
Visualización del Enrichment Map:
> gsea_go_sim <- pairwise_termsim(gsea_go)
> emapplot(gsea_go_sim, showCategory=10, color="pvalue", cluster.params=list(cluster=TRUE, legend=TRUE))
Visualización del Go plot:
> goplot(gsea_go, showCategory=10, cex=0.5)
Visualización del Category Netplot:
> cnetplot(gsea_go, color.params = list(foldChange = gsea_list), showCategory = 10)
Visualización del Ridgeplot:
> ridgeplot(gsea_go, showCategory = 10, label_format = 50, fill="p.adjust")+labs(x = "enrichment distribution")
Visualización del GSEA plot:
> # El usuario deberá seleccionar el conjunto de genes a evaluar.
> # En el ejemplo utilizamos el primer conjunto de genes.
> gseaplot(gsea_go, by = "all", title = gsea_go$Description[1], geneSetID = 1)
Para realizar el análisis GSEA con anotación en KEGG utilizaremos la función gseKEGG(). Es importante confirmar que nuestra lista de genes utiliza el identificador ENTREZID, ya que es uno de los identificadores compatibles con la función gseKEGG. En caso de que nuestra lista de genes no esté identificada con ENTREZID, se tendrá que realizar la conversión utilizando la función bitr() de clusterProfiler. Como se ha comentado anteriormente, en este trabajo la lista de genes de partida está codificada en ENTREZ ID, por lo que no se requiere de ninguna acción relacionada.
El objeto resultante almacena las vías KEGG enriquecidas, la cuenta de genes anotados en ellas y los valores estadísticos que justifican que dichas vías se encuentran significativamente enriquecidas.
> gsea_kegg <- gseKEGG(geneList = gsea_list,
+ organism = organismo_kegg,
+ minGSSize = 20,
+ pvalueCutoff = 0.05,
+ eps=0,
+ pAdjustMethod = "BH",
+ verbose= FALSE,
+ keyType = "ncbi-geneid") #Es el EntrezID
Podemos guardar los resultados obtenidos del análisis GSEA-KEGG en un archivo csv.
> gsea_kegg_results <- data.frame(gsea_kegg)
> head(gsea_kegg_results)[2:8]
Description
mmu03010 Ribosome - Mus musculus (house mouse)
mmu05168 Herpes simplex virus 1 infection - Mus musculus (house mouse)
mmu05323 Rheumatoid arthritis - Mus musculus (house mouse)
mmu04110 Cell cycle - Mus musculus (house mouse)
mmu04060 Cytokine-cytokine receptor interaction - Mus musculus (house mouse)
mmu03008 Ribosome biogenesis in eukaryotes - Mus musculus (house mouse)
setSize enrichmentScore NES pvalue p.adjust
mmu03010 133 0.6330776 2.345899 4.235348e-13 1.262134e-10
mmu05168 383 -0.3706738 -1.821980 5.158712e-10 7.686480e-08
mmu05323 66 0.6825391 2.307978 1.107689e-09 1.100304e-07
mmu04110 151 0.5573353 2.082007 2.227016e-09 1.659127e-07
mmu04060 175 0.5365291 2.042815 4.822689e-09 2.874323e-07
mmu03008 75 0.6235889 2.149150 1.470316e-07 7.302570e-06
qvalue
mmu03010 1.034317e-10
mmu05168 6.299058e-08
mmu05323 9.016974e-08
mmu04110 1.359652e-07
mmu04060 2.355503e-07
mmu03008 5.984445e-06
> write.csv(gsea_kegg_results, "Results/clusterProfiler_GSEA_KEGG.csv")
Visualización del Upset plot:
> upsetplot(gsea_kegg)
Visualización del Dotplot.
> dotplot(gsea_kegg, showCategory = 10, font.size=8, title = "Enriched Pathways" , split=".sign") + facet_grid(.~.sign)
Visualización del Enrichment Map:
> gsea_kegg_sim <- pairwise_termsim(gsea_kegg)
> emapplot(gsea_kegg_sim, showCategory=10, color="pvalue", cluster.params=list(cluster=TRUE, legend=TRUE))
Visualización del Category Netplot:
> cnetplot(gsea_kegg, color.params = list(foldChange = gsea_list), showCategory = 5)
Visualización del Ridge plot:
> ridgeplot(gsea_kegg, showCategory = 10, label_format = 50, fill="p.adjust")+labs(x = "enrichment distribution")
Visualización del GSEA plot:
> # El usuario deberá seleccionar el conjunto de genes a evaluar.
> # En el ejemplo utilizamos el primer conjunto de genes.
> gseaplot(gsea_kegg, by = "all", title = gsea_kegg$Description[1], geneSetID = 1)
Visualización del Pathview: Estas instrucciones crearán diferentes archivos con información de las vías KEGG enriquecidas. El usuario deberá seleccionar qué vía de las obtenidas en el análisis quiere visualizar.
> names(gsea_kegg@geneSets)[1]
[1] "mmu00010"
> dme <- pathview(gene.data=gsea_list, pathway.id="mmu03010", species = organismo_kegg)
> #En la aplicación se podría crear una lista de selección en base a los datos del análisis gsea_kegg@geneSets
> dme <- pathview(gene.data=gsea_list, pathway.id=gsea_kegg@result$ID[1], species = organismo_kegg)
También podríamos navegar directamente a la página web de KEGG para visualizar la vía del término seleccionado.
> browseKEGG(gsea_kegg, gsea_kegg@result$ID[1])
GSEA-KEGG Pathway para mmu00010
Para realizar el análisis GSEA con anotación con los términos Reactome utilizaremos la función gsePathway() del paquete ReactomePA. De nuevo, la función requiere que nuestra lista de genes sea un vector de identificadores ENTREZID. En caso de que nuestra lista de genes no esté identificada con ENTREZID, se tendrá que realizar la conversión utilizando la función bitr() de clusterProfiler.
Como se ha comentado anteriormente, en este trabajo la lista de genes de partida está codificada en ENTREZ ID, por lo que no se requiere de ninguna acción relacionada.
El objeto resultante almacena las vías Reactome enriquecidas, la cuenta de genes anotados en ellas y los valores estadísticos que justifican que dichas vías se encuentran significativamente enriquecidas..
> gsea_reactome <- gsePathway(gene=gsea_list,
+ organism = organismo_reactome,
+ minGSSize = 20,
+ pvalueCutoff = 0.2,
+ eps=0,
+ pAdjustMethod = "BH",
+ verbose= FALSE)
Podemos guardar los resultados obtenidos del análisis GSEA-REACTOME en un archivo csv.
> gsea_reactome_results <- data.frame(gsea_reactome)
> head(gsea_reactome_results)[2:7]
Description
R-MMU-6791226 Major pathway of rRNA processing in the nucleolus and cytosol
R-MMU-72312 rRNA processing
R-MMU-8868773 rRNA processing in the nucleus and cytosol
R-MMU-72613 Eukaryotic Translation Initiation
R-MMU-72737 Cap-dependent Translation Initiation
R-MMU-72706 GTP hydrolysis and joining of the 60S ribosomal subunit
setSize enrichmentScore NES pvalue p.adjust
R-MMU-6791226 166 0.6659914 2.531681 2.054468e-20 4.697884e-18
R-MMU-72312 166 0.6659914 2.531681 2.054468e-20 4.697884e-18
R-MMU-8868773 166 0.6659914 2.531681 2.054468e-20 4.697884e-18
R-MMU-72613 109 0.6746799 2.458890 2.406093e-14 3.301159e-12
R-MMU-72737 109 0.6746799 2.458890 2.406093e-14 3.301159e-12
R-MMU-72706 102 0.6761809 2.455942 1.228554e-13 1.404647e-11
> write.csv(gsea_reactome_results, "Results/clusterProfiler_GSEA_reactome.csv")
Visualización del Upset plot:
> upsetplot(gsea_reactome)
Visualización del Dotplot.
> dotplot(gsea_reactome, showCategory = 5, font.size=8, title = "Enriched Pathways" , split=".sign") + facet_grid(.~.sign)
Visualización del Enrichment Map:
> gsea_reactome_sim <- pairwise_termsim(gsea_reactome)
> emapplot(gsea_reactome_sim, showCategory=10, color="pvalue", cluster.params=list(cluster=TRUE, legend=TRUE))
Visualización del Category Netplot:
> cnetplot(gsea_reactome, color.params = list(foldChange = universe_list), showCategory = 5)
Visualización del Ridge plot:
> ridgeplot(gsea_reactome, showCategory = 10, label_format = 50, fill="p.adjust")+labs(x = "enrichment distribution")
Visualización del GSEA plot:
> # El usuario deberá seleccionar el conjunto de genes a evaluar.
> # En el ejemplo utilizamos el primer conjunto de genes.
> gseaplot(gsea_reactome, by = "all", title = gsea_reactome$Description[1], geneSetID = 1)
Visualización del viewPathway: Esta función permite visualizar la via Reactome indicada por el usuario.
> viewPathway(gsea_reactome@result$Description[1], organism="mouse")
El análisis de Impacto de Vías de Señalización (SPIA) es un método de análisis de enriquecimiento basado en la topología de las vías de señalización. A diferencia de los métodos anteriores, SPIA integra los datos de expresión diferencial con la topología de las vías para identificar las categorías afectadas. Es decir, el método SPIA no solo se considera los datos de expresión diferencial de la lista de genes de partida, sino también la interconexión entre los genes y su posición relativa en la vía.
El resultado del análisis SPIA es una tabla que muestra las vías significativamente desreguladas basadas en la sobre-representación y la acumulación de perturbaciones en la señalización.
Para utilizar el método SPIA partiremos de los resultados obtenidos en un estudio de expresión diferencial, en el ejemplo, nuestra topTable inicial. El análisis requiere extraer de esta matriz los valores que definen el cambio de expresión génica de los genes seleccionados (logFC) y nombrar los genes con identificador ENTREZID. El análisis evaluará la sobrerepresentación de las vías relacionadas con los genes seleccionados respecto a una lista de referencia.
La lista de genes de referencia también debe utilizar identificadores ENTREZID. Sería equivalente a la lista definida como universo en el análisis ORA. En base a la necesidad del análisis, la lista de referencia puede ser la lista de genes de la plataforma utilizada en el experimento, lo que podríamos asignar al genoma del organismo, o la propia TopTable inicial. En nuestro caso, utilizaremos la topTable inicial.
> all_spia <- as.character(df$EntrezID)
> # Eliminamos posibles valores NA
> all_spia<-na.omit(all_spia)
> # Eliminamos posibles valores duplicados
> all_spia <- all_spia[which(duplicated(all_spia) == F)]
A continuación, seleccionamos la lista de genes para analizar. En este caso, filtramos la topTable inicial en base al umbral de adj.P.Val y logFC seleccionado por el usuario.
> # De la lista de genes original, seleccionamos los de interés en base a valor adj.P.val (padj < 0.05)
> # En la aplicación se añadirán selectores para filtrar los datos
> sig_genes_spia = subset(df, adj.P.Val < 0.05)
> # De la lista anterior, extraemos los valores de logFC para filtrar posteriormente
> genes_spia <- sig_genes_spia$logFC
> # Nombramos el vector anterior con los identificadores
> names(genes_spia) <- sig_genes_spia$EntrezID
> # omitimos NA values
> genes_spia <- na.omit(genes_spia)
> # Eliminamos posibles valores duplicados
> genes_spia <- genes_spia[which(duplicated(names(genes_spia)) == F)]
> # Finalmente filtramos por el valor de logFC deseado. En este caso nos centraremos en los genes up/down-regulated
> # abs(logFC) > 2
> genes_spia <- genes_spia[abs(genes_spia) > 2]
Realizamos el análisis con la función spia() del paquete SPIA. Se debe indicar en el argumento de el vector de los genes en evaluación con los resultados de expresión génica (logFC) y en argumento all la lista de genes de referencia. Como se ha comentado anteriormente, esta función sólo admite la identificción en ENTREZID.
Uno de los parámetros configurables en el análisis es el número de iteraciones bootstrap nB. Suele ser mayor a 100, recomendando un valor de 2000. Con el objetivo de reducir el tiempo de procesado, en el ejemplo se utilizará un valor de 100.
El objeto resultante almacena las categorías KEGG enriquecidas en el subconjunto de genes. La tabla de resultados muestra la siguiente información para cada vía:
pSize: el número de genes en la vía.
NDE: el número de genes diferencialmente expresados en cada categoría.
tA: la acumulación total observada de perturbaciones en la vía.
pNDE: la probabilidad de observar al menos NDE genes en la vía utilizando un modelo hipergeométrico (similar a ORA).
pPERT: la probabilidad de observar una acumulación total de perturbaciones igual a tA de manera aleatoria.
pG: el valor p obtenido combinando pNDE y pPERT.
pGFdr y pGFWER son los valores p globales ajustados FDR y Bonferroni, respectivamente.
Status: indica la dirección en la que se perturba la vía (activada o inhibida).
KEGGLINK proporciona un enlace web al sitio web de KEGG que muestra la ruta de la vía.
> spia_kegg <- spia(de=genes_spia, all=all_spia, nB=100, organism=organismo_kegg, verbose = FALSE)
Podemos guardar los resultados obtenidos del análisis SPIA en un archivo csv.
> spia_kegg_results <- data.frame(spia_kegg)
> head(spia_kegg_results)[,1:6]
Name ID pSize NDE pNDE tA
1 Staphylococcus aureus infection 05150 30 8 5.127886e-07 36.13560
2 Systemic lupus erythematosus 05322 67 10 5.910330e-06 21.50081
3 Complement and coagulation cascades 04610 43 7 8.561225e-05 43.83694
4 Osteoclast differentiation 04380 101 6 4.100137e-02 23.47751
5 ECM-receptor interaction 04512 83 4 1.538758e-01 -12.48262
6 Pertussis 05133 61 7 7.884269e-04 14.81044
> head(spia_kegg_results)[,7:11]
pPERT pG pGFdr pGFWER Status
1 1e-04 1.266267e-09 1.392894e-07 1.392894e-07 Activated
2 1e-04 1.314998e-08 7.232489e-07 1.446498e-06 Activated
3 1e-04 1.675947e-07 6.145140e-06 1.843542e-05 Activated
4 1e-04 5.496025e-05 1.511407e-03 6.045627e-03 Activated
5 1e-04 1.859119e-04 4.090062e-03 2.045031e-02 Inhibited
6 6e-02 5.184166e-04 9.504304e-03 5.702583e-02 Activated
> head(spia_kegg_results)[,12]
[1] "http://www.genome.jp/dbget-bin/show_pathway?mmu05150+12266+12259+12260+12262+12268+14960+14961+19204"
[2] "http://www.genome.jp/dbget-bin/show_pathway?mmu05322+12268+12266+12259+12260+12262+21926+14960+14961+319189+11472"
[3] "http://www.genome.jp/dbget-bin/show_pathway?mmu04610+14061+18787+12266+12259+12260+12262+12268"
[4] "http://www.genome.jp/dbget-bin/show_pathway?mmu04380+12978+12977+22177+16822+21926+17970"
[5] "http://www.genome.jp/dbget-bin/show_pathway?mmu04512+12833+12834+21827+15366"
[6] "http://www.genome.jp/dbget-bin/show_pathway?mmu05133+20311+12259+12260+12262+12268+12266+21926"
> write.csv(spia_kegg_results, "Results/SPIA.csv")
Visualización del PlotP.
Se debe definir el parámetro threshold, un valor numérico entre 0 y 1 utilizado como umbral de significación de la categoría. Normalmente se utiliza un valor de 0.05.
> #Utilizando la función integrada en SPIA plotP obtenemos el siguiente error, que parece ser un fallo de programación de la función:
> #Error in if (sum(oky) > 0) { : missing value where TRUE/FALSE needed
>
> #Con el objetivo de resolver el error se modifica ligeramente el código base de la función plotP() encontada en https://rdrr.io/bioc/SPIA/src/R/plotP.R
>
> getP2<-function(pG,combine="fisher"){
+ #given a pG returns two equal p-values such as combfunc(p1,p2)=pG
+ if(combine=="fisher"){
+ ch=qchisq(pG,4,lower.tail = FALSE)
+ return(sqrt(exp(-ch/2)))
+ }
+
+ if(combine=="norminv"){
+ return(pnorm(qnorm(pG)*sqrt(2)/2))
+ }
+ }
>
> plotP2<-function(x,threshold=0.05){
+
+ if(class(x)!="data.frame" | dim(x)[1]<1 | !all(c("ID","pNDE","pPERT","pG","pGFdr","pGFWER")%in%names(x)))
+ {
+ stop("plotP can be applied only to a dataframe produced by spia function!!!")
+ }
+
+
+ if(threshold<x[1,"pGFdr"]){
+ msg<-paste("The threshold value should be",x[1,"pGFdr"],"or higher!!!");
+ stop(msg);
+ }
+
+ pb<-x[,"pPERT"]
+ ph<-x[,"pNDE"]
+
+ #determine what combine method was used to convert ph and pb into pG
+ combinemethod=ifelse(sum(combfunc(pb,ph,"fisher")==x$pG)>sum(combfunc(pb,ph,"norminv")==x$pG),"fisher","norminv")
+
+
+ okx<-(ph<1e-6)
+ oky<-(pb<1e-6)
+
+ ph[ph<1e-6]<-1e-6
+ pb[pb<1e-6]<-1e-6
+
+ plot(-log(ph),-log(pb),xlim=c(0,max(c(-log(ph),-log(pb))+1,na.rm=TRUE)),
+ ylim=c(0,max(c(-log(ph),-log(pb)+1),na.rm=TRUE)),pch=19,main="SPIA two-way evidence plot",cex=1.5,
+ xlab="-log(P NDE)",ylab="-log(P PERT)")
+ tr<-threshold/dim(na.omit(x))[1]
+ abline(v=-log(tr),lwd=1,col="red",lty=2)
+ abline(h=-log(tr),lwd=1,col="red",lty=2)
+
+ if(combinemethod=="fisher"){
+ points(c(0,-log(getP2(tr,"fisher")^2)),c(-log(getP2(tr,"fisher")^2),0),col="red",lwd=2,cex=0.7,type="l")
+ }else{
+ somep1=exp(seq(from=min(log(ph)),to=max(log(ph)),length=200))
+ somep2=pnorm(qnorm(tr)*sqrt(2)-qnorm(somep1))
+ points(-log(somep1),-log(somep2),col="red",lwd=2,cex=0.7,type="l")
+ }
+
+ oks<-x[,"pGFWER"]<=threshold
+ trold=tr
+ tr<-max(x[,"pG"][x[,"pGFdr"]<=threshold])
+ if(tr<=trold){tr=trold*1.03}
+
+ if(combinemethod=="fisher"){
+ points(c(0,-log(getP2(tr,"fisher")^2)),c(-log(getP2(tr,"fisher")^2),0),col="blue",lwd=2,cex=0.7,type="l")
+ }else{
+ somep1=exp(seq(from=min(log(ph)),to=max(log(ph)),length=200))
+ somep2=pnorm(qnorm(tr)*sqrt(2)-qnorm(somep1))
+ points(-log(somep1),-log(somep2),col="blue",lwd=2,cex=0.7,type="l")
+ }
+
+ abline(v=-log(tr),lwd=1,col="blue",lty=2)
+ abline(h=-log(tr),lwd=1,col="blue",lty=2)
+ text(-log(ph)[oks]+0.70,-log(pb)[oks],labels=as.vector(x$ID)[oks],cex=0.65)
+ oks2<-x[,"pGFdr"]<=threshold
+ points(-log(ph)[oks2],-log(pb)[oks2],pch=19,col="blue",cex=1.5)
+ points(-log(ph)[oks],-log(pb)[oks],pch=19,col="red",cex=1.5)
+
+ text(-log(ph)[oks2]+0.70,-log(pb)[oks2],labels=as.vector(x$ID)[oks2],cex=0.65)
+
+ if (sum(okx) > 0) {
+ points(-log(ph)[okx] - 0.12, -log(pb)[okx], pch = "|", col = "black", cex = 1.5)
+ }
+ if (sum(!is.na(oky) & oky > 0)) { #Se utiliza la función is.na() para verificar si oky tiene valores NA
+ points(-log(ph)[oky[!is.na(oky) & oky > 0]], -log(pb)[oky[!is.na(oky) & oky > 0]] - 0.12, pch = "_", col = "black", cex = 1.5)
+ }
+
+ }
>
> plotP2(spia_kegg, threshold=0.05)
Visualización las categorías KEGG. Podemos navegar directamente a la página web de KEGG para visualizar la vía del término seleccionado.
> pathway_selected<-1 #A seleccionar por el usuario. En este caso seleccionamos la primera.
> browseURL(spia_kegg[pathway_selected, 12])
SPIA-KEGG Pathway para 05150
El código presentado en este trabajo proviene de la revisión, evaluación e implementación de diferentes fuentes bibliográficas:
Ten Years of Pathway Analysis: Current Approaches and Outstanding Challenges
Statistical analysis and visualization of functional profiles for genes and gene clusters
Seamless navigation through combined results of set- & network-based enrichment analysis
Biomedical Knowledge Mining using GOSemSim and clusterProfiler
Finalmente, se reportan las versiones de los paquetes utilizados en el informe:
> packageVersion("clusterProfiler")
[1] '4.6.2'
> packageVersion("enrichplot")
[1] '1.18.3'
> packageVersion("ReactomePA")
[1] '1.42.0'
> packageVersion("ggplot2")
[1] '3.4.1'
> packageVersion("DOSE")
[1] '3.24.2'
> packageVersion("pathview")
[1] '1.38.0'
> packageVersion("SPIA")
[1] '2.50.0'
> packageVersion(organismo_go)
[1] '3.16.0'